/*
  File:    AdvectTumor.cpp
  Purpose: These routines construct normal velocities and advect the tumor. 

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <iostream>

using namespace std;

#include "../Matrix.h"
#include "Gradients.h"
#include "../ImprovedMath.h"
#include "GaussianFilter.h"
#include "GeometricData.h"
#include "Curvature.h"
#include "WENO.h"

#include "AdvectTumor.h"

#include "Polynomial.h"
#include "../ParseText.h"

Polynomial ECMtaxisPolynomial;
bool ECMtaxisPolynomialSet = false;
double ECMtaxisLeftCutoff;
double ECMtaxisRightCutoff;
double MinECMtaxis;
double MaxECMtaxis;

bool CreateECMtaxisPolynomial( void )
{
 if( ECMtaxisPolynomialSet ) 
 { return false; }

 extern char* ParameterFile;

 ECMtaxisLeftCutoff = RetrieveDoubleParameter(ParameterFile,"ECMtaxisLeftCutoff");
 ECMtaxisRightCutoff = RetrieveDoubleParameter(ParameterFile,"ECMtaxisRightCutoff");
 MinECMtaxis = RetrieveDoubleParameter(ParameterFile,"MinECMtaxis");
 MaxECMtaxis = RetrieveDoubleParameter(ParameterFile,"MaxECMtaxis");

 ECMtaxisPolynomial.SetDegree(4);
 for( int i=0 ; i <= ECMtaxisPolynomial.GetDegree() ; i++ )
 { ECMtaxisPolynomial.SetCoefficient( i , 1.0 ); }
 
 Polynomial ChiPrime = Derivative( ECMtaxisPolynomial );
 
 Matrix Coeff = zeros( ECMtaxisPolynomial.GetDegree()+1,ECMtaxisPolynomial.GetDegree()+1 );
 Matrix RHS = zeros(1,ECMtaxisPolynomial.GetDegree()+1 );
 
 for( int j=0; j <= ECMtaxisPolynomial.GetDegree() ; j++ )
 { *Coeff(0,j) = ECMtaxisPolynomial.GetCoefficient(j)*IntPow( ECMtaxisLeftCutoff , j ); }
 *RHS(0) = MinECMtaxis;
 *Coeff(1,0) = 0.0;
 for( int j=1; j <= ECMtaxisPolynomial.GetDegree() ; j++ )
 { *Coeff(1,j) = ChiPrime.GetCoefficient(j-1)*IntPow( ECMtaxisLeftCutoff , j-1 ); }
 *RHS(1) = 0.0;
 
 double ECMchiMiddle = 0.5*(ECMtaxisLeftCutoff + ECMtaxisRightCutoff);
 
 for( int j=0; j <= ECMtaxisPolynomial.GetDegree() ; j++ )
 { *Coeff(2,j) = ECMtaxisPolynomial.GetCoefficient(j)*IntPow( ECMchiMiddle , j ); }
 *RHS(2) = MaxECMtaxis;
 
 for( int j=0; j <= ECMtaxisPolynomial.GetDegree() ; j++ )
 { *Coeff(3,j) = ECMtaxisPolynomial.GetCoefficient(j)*IntPow( ECMtaxisRightCutoff , j ); }
 *RHS(3) = MinECMtaxis;
 *Coeff(4,0) = 0.0;
 for( int j=1; j <= ECMtaxisPolynomial.GetDegree() ; j++ )
 { *Coeff(4,j) = ChiPrime.GetCoefficient(j-1)*IntPow( ECMtaxisRightCutoff , j-1 ); }
 *RHS(4) = 0.0; 
 
 Matrix FinalChiCoefficients = GaussianEliminationPartialPivoting( Coeff, RHS  );
 
 for( int i=0; i <= ECMtaxisPolynomial.GetDegree() ; i++ )
 { ECMtaxisPolynomial.SetCoefficient( i , *FinalChiCoefficients(i) ); }
 ChiPrime = Derivative( ECMtaxisPolynomial );
 
 ECMtaxisPolynomialSet = true;
 return true;
}

double ECMtaxisFunction( double ecm )
{
 if( ecm < ECMtaxisLeftCutoff || ecm > ECMtaxisRightCutoff ) 
 { return MinECMtaxis; }
 return ECMtaxisPolynomial( ecm ); 
}

bool UpdateECMtaxisCoefficient( void )
{
 CreateECMtaxisPolynomial();

 extern Matrix ECMtaxisCoefficient;
 extern int Xnodes, Ynodes;
 extern double dx;
 extern Matrix X,Y;
 extern double _EPS;
 
 ECMtaxisCoefficient = ones( Xnodes, Ynodes );

 extern Matrix Phi, PhiNecrotic, ECM;
 extern double FunctionSmoothingLevel;
 Matrix ECMsmoothed = GaussianFilter(X,Y,ECM,FunctionSmoothingLevel*dx);
 
 for( int i=0 ; i < Xnodes ; i++ )
 {
  for( int j=0 ; j < Ynodes ; j++ )
  {
   *ECMtaxisCoefficient(i,j) = ECMtaxisFunction( *ECMsmoothed(i,j) );
  }
 }
/* 
 extern double FunctionSmoothingLevel;

 ECMtaxisCoefficient = GaussianFilter(X,Y,ECMtaxisCoefficient,FunctionSmoothingLevel*dx);
*/
 return true;
}

bool UpdateInteriorECMtaxisCoefficient( void )
{
 CreateECMtaxisPolynomial();

 extern Matrix ECMtaxisCoefficient;
 extern int Xnodes, Ynodes;
 extern double dx;
 extern Matrix X,Y;
 extern double _EPS;
 
 ECMtaxisCoefficient = ones( Xnodes, Ynodes );
 
 extern Matrix Phi, PhiNecrotic, ECM;
 extern double FunctionSmoothingLevel;
 Matrix ECMsmoothed = GaussianFilter(X,Y,ECM,FunctionSmoothingLevel*dx); 
 
 for( int i=0 ; i < Xnodes ; i++ )
 {
  for( int j=0 ; j < Ynodes ; j++ )
  {
   if( *Phi(i,j) > _EPS || *PhiNecrotic(i,j) <= _EPS )
   { *ECMtaxisCoefficient(i,j) = 0.0; }
   else
   { *ECMtaxisCoefficient(i,j) = ECMtaxisFunction( *ECMsmoothed(i,j) ); }
  }
 }
 
/* 
 extern double FunctionSmoothingLevel;

 ECMtaxisCoefficient = GaussianFilter(X,Y,ECMtaxisCoefficient,FunctionSmoothingLevel*dx);
*/
 return true;
}


double BilinearInterpolation( Matrix &Xin, Matrix &Yin, Matrix& Data, double x, double y )
{
 int i=0;
 while( *Xin(i) < x && i < Xin.TellCols()-2 )
 { i++; }
 i--;
 int j=0;
 while( *Yin(j) < y && j < Yin.TellCols()-2 )
 { j++; }
 j--;
 
 if( x < *Xin(0) || x > *Xin( Xin.TellCols()-2 ) )
 {
  cout << "BilinearInterpolation error: requested out-of-range point!" << endl;
  return 0.0; 
 }
 if( y < *Yin(0) || y > *Yin( Yin.TellCols()-2 ) )
 {
  cout << "BilinearInterpolation error: requested out-of-range point!" << endl;
  return 0.0;
 }
 
 double dXinterpolationBottom = ( *Xin(i+1) - x )*( *Data(i,j)   ) 
                              + ( x - *Xin(i)   )*( *Data(i+1,j) );
 dXinterpolationBottom /= ( *Xin(1)-*Xin(0) );     

 double dXinterpolationTop = ( *Xin(i+1) - x )*( *Data(i,j+1)   ) 
                           + ( x - *Xin(i)   )*( *Data(i+1,j+1) );
 dXinterpolationTop /= (*Xin(1)-*Xin(0));   
     
 double dOut = ( *Yin(j+1) - y )*( dXinterpolationBottom  )
             + ( y - *Yin(j)   )*( dXinterpolationTop     );
 dOut /= (*Yin(1)-*Yin(0));
 return dOut;
}

Matrix BilinearExtension( Matrix& Input )
{ 
 extern Matrix Phi;
 extern Matrix X;
 extern Matrix Y;

 Matrix Output = zeros( Input.TellRows() , Input.TellCols() );
  
 int i,j;

 int Count = 0;
 int TotalCount = 0;
 
 for( j=1 ; j < Output.TellCols()-1; j++ )
 {
  for( i=1 ; i < Output.TellRows()-1 ; i++ )
  {
   if( *Phi(i,j) <= 2.0 )
   {
    GeometricDataPacket GDP = AdaptiveGeometricData(X,Y,Phi,*X(i),*Y(j));
    
    double nx = GDP.NormalX; 
    double ny = GDP.NormalY; 
    double norm = sqrt( square(nx) +square(ny) );
    
    double TargetX = *X(i) - *Phi(i,j) * nx / norm;
    double TargetY = *Y(j) - *Phi(i,j) * ny / norm;
   
    double phi = BilinearInterpolation(X,Y,Phi,TargetX,TargetY);
    if( fabs( phi ) > square(dx) )
    {
     Count++;
    }
    
    *Output(i,j) = BilinearInterpolation(X,Y,Input,TargetX,TargetY);
    TotalCount++;
   }
  }
 }
 cout << "bad at " << (double) Count * 100.0 /(double) TotalCount << "% of nodes" << endl;
 
 return Output;
}

Matrix BilinearExtension( Matrix& LevelSet, Matrix& Input )
{ 
 extern Matrix X;
 extern Matrix Y;

 Matrix Output = zeros( Input.TellRows() , Input.TellCols() );
  
 int i,j;

 for( j=1 ; j < Output.TellCols()-1; j++ )
 {
  for( i=1 ; i < Output.TellRows()-1 ; i++ )
  {
   if( *LevelSet(i,j) <= 2.0 )
   {
    GeometricDataPacket GDP = AdaptiveGeometricData(X,Y,LevelSet,*X(i),*Y(j));
    double TargetX = *X(i) - *LevelSet(i,j) * GDP.NormalX;
    double TargetY = *Y(j) - *LevelSet(i,j) * GDP.NormalY;
    *Output(i,j) = BilinearInterpolation(X,Y,Input,TargetX,TargetY);
   }
  }
 }
 return Output;
}


bool ConstructVelocities( void )
{
 extern Matrix PressureGradientX;
 extern Matrix PressureGradientY;
 extern Matrix Pressure;
 extern Matrix Speed;
 extern Matrix NecroticSpeed;
 extern Matrix Mobility;
 extern Matrix X;
 extern Matrix Y;
 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 extern int Xnodes;
 extern int Ynodes;
 extern double dx;

 PressureGradientX = zeros( Xnodes, Ynodes );
 PressureGradientY = zeros( Xnodes, Ynodes );

 double FilterStandardDeviation = 5.0*dx;

 bool Output = CalculateGradients( Pressure, &PressureGradientX, &PressureGradientY );
 
 extern Matrix ECM;
 extern Matrix ECMgradientX;
 extern Matrix ECMgradientY;
 
 ECMgradientX = zeros(Xnodes,Ynodes);
 ECMgradientY = zeros(Xnodes,Ynodes);
 
 extern double FunctionSmoothingLevel;
 Matrix SmoothedECM = GaussianFilter(X,Y,ECM, FunctionSmoothingLevel*dx);
 
 extern Matrix ECMtaxisCoefficient;
 Output &= UpdateInteriorECMtaxisCoefficient();
 
 Output &= CalculateGradients( SmoothedECM, &ECMgradientX, &ECMgradientY );
 
 Speed = zeros( Xnodes, Ynodes );
 

 int i,j; 
 for( j=1 ; j < Ynodes-1 ; j++ )
 {
  for( i=1 ; i < Xnodes-1 ; i++ )
  {
   if( *Phi(i,j) <= 2.0 )
   {   
    GeometricDataPacket GDP = AdaptiveGeometricData(X,Y,Phi,*X(i),*Y(j));
    *Speed(i,j) = -*Mobility(i,j)*(   *PressureGradientX(i,j) * GDP.NormalX 
                                    + *PressureGradientY(i,j) * GDP.NormalY );
    *Speed(i,j) += *ECMtaxisCoefficient(i,j)*(   *ECMgradientX(i,j) * GDP.NormalX 
                                               + *ECMgradientY(i,j) * GDP.NormalY );	
   }                                    
  }
 }

 ExtrapolateInteriorSpeed( Phi, &Speed );

 Speed = BilinearExtension( Speed );
 Speed = GaussianFilter(X,Y,Speed, FilterStandardDeviation); // 
 Speed = BilinearExtension( Speed );
 
 // if there is an error, it must occur somewhere around here. 1-2-2007

 // redo the geometric information to suit the necrotic level set function, if there is a necrotic core
 
 bool Exists = false;
 for( j=1 ; j < Ynodes-1 ; j++ )
 {
  for( i=1 ; i < Xnodes-1 ; i++ )
  {
   if( *PhiNecrotic(i,j) <= _EPS )
   { Exists = true; }
  }
 }

 NecroticSpeed = zeros(Xnodes,Ynodes);
 if( Exists )
 {
  GenerateGeometricGlobals( X,Y,PhiNecrotic );
  for( j=1 ; j < Ynodes-1 ; j++ )
  {
   for( i=1 ; i < Xnodes-1 ; i++ )
   {
    if( *PhiNecrotic(i,j) <= 2.0 )
    {   
     GeometricDataPacket GDP = AdaptiveGeometricData(X,Y,PhiNecrotic,*X(i),*Y(j));
     *NecroticSpeed(i,j) = -*Mobility(i,j)*(   *PressureGradientX(i,j) * GDP.NormalX 
                                             + *PressureGradientY(i,j) * GDP.NormalY );
     *NecroticSpeed(i,j) += *ECMtaxisCoefficient(i,j)*(  *ECMgradientX(i,j) * GDP.NormalX 
                                                        +*ECMgradientY(i,j) * GDP.NormalY );	
											 
    }                                             
   }
  }
  NecroticSpeed = BilinearExtension( PhiNecrotic, NecroticSpeed );
  NecroticSpeed = GaussianFilter(X,Y,NecroticSpeed,FilterStandardDeviation);
  NecroticSpeed = BilinearExtension( PhiNecrotic, NecroticSpeed );
  
  // don't forget to switch back to the tumor normals and geometric data!
  GenerateGeometricGlobals( X,Y,Phi );
 }

 return Output;
}

double AdvectTumor( void )
{
 extern Matrix Speed;
 extern Matrix Phi;
 extern int Xnodes;
 extern int Ynodes;
 extern Matrix PhiNecrotic;
 extern Matrix NecroticSpeed;
 extern double t;
 extern double NextOutputTime;

 Matrix PhiOperator = WENO5(Speed, Phi);
 Matrix PhiNecroticOperator = zeros(Xnodes,Ynodes);
 
 bool Exists = false;
 int i,j;
 for( j=1 ; j < Ynodes-1 ; j++ )
 {
  for( i=1 ; i < Xnodes-1 ; i++ )
  {
   if( *PhiNecrotic(i,j) <= _EPS )
   { Exists = true; }
  }
 }

 if( Exists )
 {
  PhiNecroticOperator = WENO5(NecroticSpeed, PhiNecrotic);
 }

 double MaxSpeedA = Speed.InfinityNorm();
 double MaxSpeedB = NecroticSpeed.InfinityNorm();
 double MaxSpeed = MAX( MaxSpeedA, MaxSpeedB );
 
 double DeltaT = 0.5 * dx / MaxSpeed;
//double DeltaT=0.05;//ADDED TEMP
 
 // new on 11-16-2006
 if( DeltaT > 0.05 )
 { DeltaT = 0.05; }

 if( t + DeltaT > NextOutputTime )
 { DeltaT = NextOutputTime - t; }
  
 Phi = Phi - DeltaT*PhiOperator;
 if( Exists )
 { PhiNecrotic = PhiNecrotic - DeltaT*PhiNecroticOperator; }
 

 return DeltaT; 
}

